1. Introducción al estudio

Hemos escogido un dataset que recoge ciertas características de unas casas en el área del condado de King, del estado de Whasington junto al valor por el que se vendieron. En el conjunto tenemos un total de 21613 observaciones, es decir viviendas.

Nuestro objetivo será usar esta muestra para poder predecir el valor de una casa en este área en funcion de las variables estudiadas.

Las variables que nos ocupan y que se presentan en el dataset son:

  • id – Identificador único de cada casa
  • date – Fecha en la que se vendió la casa
  • price – Precio
  • bedrooms – Número de habitaciones
  • bathrooms – Número de baños
  • sqft_living – Pies cuadrados útiles de la casa
  • sqft_lot – Pies cuadrados del terreno en el que se sitúa la casa
  • floors – Número de plantas
  • waterfront – Si el apartamento tiene o no vistas al mar
  • view – Cómo de buena son las vistas desde la propiedad
  • condition – Condición de la casa
  • grade – Nivel de construcción y diseño de la casa
  • sqft_above – Pies cuadrados del interior de la casa sin contar el sótano
  • sqft_basement – Pies cuadrados del sótano
  • yr_built – Año en que fue construida la casa inicialmente
  • yr_renovated – Año en el que la casa fue renovada
  • zipcode – Cógido referente a la zona donde se encuentra la casa
  • lat - Latitud
  • long - Longitud
  • sqft_living15 – media de los pies cuadrados útiles de las 15 casas más cercanas
  • sqft_lot15 – media de los pies cuadrados del terreno de las 15 casas más cercanas

2. Tranformaciones en el tipo de variables

str(datos)
## 'data.frame':    21613 obs. of  21 variables:
##  $ id           : num  7.13e+09 6.41e+09 5.63e+09 2.49e+09 1.95e+09 ...
##  $ date         : chr  "20141013T000000" "20141209T000000" "20150225T000000" "20141209T000000" ...
##  $ price        : num  221900 538000 180000 604000 510000 ...
##  $ bedrooms     : int  3 3 2 4 3 4 3 3 3 3 ...
##  $ bathrooms    : num  1 2.25 NA 3 2 4.5 NA 1.5 NA 2.5 ...
##  $ sqft_living  : int  1180 2570 770 1960 1680 5420 1715 1060 1780 1890 ...
##  $ sqft_lot     : int  5650 7242 10000 5000 8080 101930 6819 9711 7470 6560 ...
##  $ floors       : num  1 2 NA 1 1 1 NA 1 NA 2 ...
##  $ waterfront   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ view         : int  0 0 NA 0 0 0 NA 0 NA 0 ...
##  $ condition    : int  3 3 3 5 3 3 3 3 3 3 ...
##  $ grade        : int  7 7 6 7 8 11 7 7 7 7 ...
##  $ sqft_above   : int  1180 2170 770 1050 1680 3890 1715 1060 1050 1890 ...
##  $ sqft_basement: int  0 400 0 910 0 1530 0 0 730 0 ...
##  $ yr_built     : int  1955 1951 1933 1965 1987 2001 1995 1963 1960 2003 ...
##  $ yr_renovated : int  0 1991 0 0 0 0 0 0 0 0 ...
##  $ zipcode      : int  98178 98125 NA 98136 98074 98053 NA 98198 NA 98038 ...
##  $ lat          : num  47.5 47.7 NA 47.5 47.6 ...
##  $ long         : num  -122 -122 NA -122 -122 ...
##  $ sqft_living15: int  1340 1690 2720 1360 1800 4760 2238 1650 1780 2390 ...
##  $ sqft_lot15   : int  5650 7639 8062 5000 7503 101930 6819 9711 8113 7570 ...

Podemos observar que hay variables que no están codificadas correctamente:

  • La variable “id” aparece como numérica,“num”, cuando realmente debería ser de tipo carácter,“chr”.
  • La variable “date” aparece como “chr”, cuando debería ser “Date”.
  • Las variables “waterfront” y “zipcode” deberían ser de tipo “factor” en vez de enteros “int”.

Por lo tanto, para modificar los tipos de variables hacemos lo siguiente:

datos$id=as.character(datos$id)
datos$date = as.Date(datos$date, "%Y%m%dT%H%M%S")
#datee=datos$date #lo usaremos más adelante, en el gráfico para estudiar la variable date
DATOS<- read.csv("kc_house_data.csv")
datee= as.Date(DATOS$date, "%Y%m%dT%H%M%S")
datos$waterfront = as.factor(datos$waterfront)
datos$zipcode = as.factor(datos$zipcode)
datos$condition=as.factor(datos$condition)
datos$date=as.factor(datos$date)
datos$view=as.factor(datos$view)

str(datos)
## 'data.frame':    21613 obs. of  21 variables:
##  $ id           : chr  "7129300520" "6414100192" "5631500400" "2487200875" ...
##  $ date         : Factor w/ 372 levels "2014-05-02","2014-05-03",..: 165 221 291 221 284 11 57 252 340 306 ...
##  $ price        : num  221900 538000 180000 604000 510000 ...
##  $ bedrooms     : int  3 3 2 4 3 4 3 3 3 3 ...
##  $ bathrooms    : num  1 2.25 NA 3 2 4.5 NA 1.5 NA 2.5 ...
##  $ sqft_living  : int  1180 2570 770 1960 1680 5420 1715 1060 1780 1890 ...
##  $ sqft_lot     : int  5650 7242 10000 5000 8080 101930 6819 9711 7470 6560 ...
##  $ floors       : num  1 2 NA 1 1 1 NA 1 NA 2 ...
##  $ waterfront   : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ view         : Factor w/ 5 levels "0","1","2","3",..: 1 1 NA 1 1 1 NA 1 NA 1 ...
##  $ condition    : Factor w/ 5 levels "1","2","3","4",..: 3 3 3 5 3 3 3 3 3 3 ...
##  $ grade        : int  7 7 6 7 8 11 7 7 7 7 ...
##  $ sqft_above   : int  1180 2170 770 1050 1680 3890 1715 1060 1050 1890 ...
##  $ sqft_basement: int  0 400 0 910 0 1530 0 0 730 0 ...
##  $ yr_built     : int  1955 1951 1933 1965 1987 2001 1995 1963 1960 2003 ...
##  $ yr_renovated : int  0 1991 0 0 0 0 0 0 0 0 ...
##  $ zipcode      : Factor w/ 70 levels "98001","98002",..: 67 56 NA 59 38 30 NA 69 NA 24 ...
##  $ lat          : num  47.5 47.7 NA 47.5 47.6 ...
##  $ long         : num  -122 -122 NA -122 -122 ...
##  $ sqft_living15: int  1340 1690 2720 1360 1800 4760 2238 1650 1780 2390 ...
##  $ sqft_lot15   : int  5650 7639 8062 5000 7503 101930 6819 9711 8113 7570 ...

Con “str(datos)” podemos ver cómo ya si tenemos las variables bien codificadas.

3. División de los datos es training y testing

Para poder crear nuestro modelo y ver su capacidad de predicción dividiremos el conjunto de datos en dos con un porcentaje de 70% y 30% para el grupo de training y test, respectivamente.

set.seed(12345)
inTraining <- createDataPartition(datos$price,
                                  p = .7, list = FALSE, times = 1)
d_training <- slice(datos, inTraining)
d_testing <- slice(datos, -inTraining)
ntrain=nrow(d_training)

4. Análisis exploratorio de datos faltantes: VIM

Veamos cuántas observaciones con datos faltantes tenemos en nuestro dataset:

ntrain-sum(complete.cases(d_training))
## [1] 68

Obtenemos que hay 68 observaciones de un total de 15131, que tienen datos faltantes. Esto se corresponde con una proporcion de 0.0043 de casos con algún dato faltante. Lo cual al ser menor de 0.03 no tenemos un problema grave. Entre las posibles técnicas a aplicar para la imputación de estos datos faltantes podría ser el asignar la mediana del resto de sus valores conocidos en el caso de variables cualitativas y en el caso de las variables cualitativas asignar la categoría más frecuente de sus conocidas.
Otra opción, debido al bajo número de observaciones con datos faltantes, sería descartar directamente los casos con NA.

Realizaremos unos gráficos para analizar posibles relaciones y estructuras de los valores faltantes.

aggr_plot=aggr(d_training,,col=c("green","red"),numbers=TRUE, sortVars=TRUE, labels=names(d_training),cex.axis=0.7,gap=3,ylab=c("Hist","Pat")) 

## 
##  Variables sorted by number of missings: 
##       Variable        Count
##           view 0.0040314586
##        zipcode 0.0022470425
##            lat 0.0022470425
##          grade 0.0011235212
##     waterfront 0.0009913423
##     sqft_above 0.0007930738
##             id 0.0006608949
##           date 0.0006608949
##    sqft_living 0.0003304474
##  sqft_living15 0.0003304474
##      bathrooms 0.0001982685
##         floors 0.0001982685
##           long 0.0001982685
##          price 0.0000000000
##       bedrooms 0.0000000000
##       sqft_lot 0.0000000000
##      condition 0.0000000000
##  sqft_basement 0.0000000000
##       yr_built 0.0000000000
##   yr_renovated 0.0000000000
##     sqft_lot15 0.0000000000

Los gráficos anteriores nos permiten descubrir rápidamente qué variables de nuestro dataset tienen mayor cantidad de datos faltantes (gráfico a la izquierda) y si puede existir algún patrón de co-ocurrencia en los datos faltantes de varias variables (gráfico a la derecha).

summary(aggr_plot)
## 
##  Missings per variable: 
##       Variable Count
##             id    10
##           date    10
##          price     0
##       bedrooms     0
##      bathrooms     3
##    sqft_living     5
##       sqft_lot     0
##         floors     3
##     waterfront    15
##           view    61
##      condition     0
##          grade    17
##     sqft_above    12
##  sqft_basement     0
##       yr_built     0
##   yr_renovated     0
##        zipcode    34
##            lat    34
##           long     3
##  sqft_living15     5
##     sqft_lot15     0
## 
##  Missings in combinations of variables: 
##                               Combinations Count     Percent
##  0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0 15063 99.55059150
##  0:0:0:0:0:0:0:0:0:0:0:0:1:0:0:0:0:0:0:0:0     2  0.01321790
##  0:0:0:0:0:0:0:0:0:1:0:0:0:0:0:0:1:1:0:0:0    31  0.20487740
##  0:0:0:0:0:0:0:0:0:1:0:1:0:0:0:0:0:0:0:0:0    17  0.11235212
##  0:0:0:0:0:1:0:0:1:0:0:0:0:0:0:0:0:0:0:1:0     5  0.03304474
##  0:0:0:0:1:0:0:1:0:1:0:0:0:0:0:0:1:1:1:0:0     3  0.01982685
##  1:1:0:0:0:0:0:0:1:1:0:0:1:0:0:0:0:0:0:0:0    10  0.06608949

Aquí podemos observar el porcentaje de valores perdidos para cada variable. Dichos porcentajes nos permiten ver que no existe ningun patrón de cocurrencia en los datos faltantes de las variables. Esto es, que no tenemos valores de variables que tienden a suceder conjuntamente.

Debido al bajo número de observaciones con datos faltantes y a que no existe ningún patrón de cocurrencia en los datos faltantes procedemos a eliminarlos de nuestro estudio.

d_training=d_training[complete.cases(d_training)==TRUE,] #eliminamos los datos faltantes
dim(d_training)
## [1] 15063    21

Ahora trabajaremos con un conjunto de datos con 15063 observaciones.

5. Estudio de las variables y sus correspondientes transformaciones

Price

El histograma es un gráfico que representa de forma bastante precisa la distribución del conjunto de datos pudiéndose observar su dispersión y/o asimetría. Por ello, realizaremos un histograma del precio para ver cómose comporta esta variable.

d_training %>% ggplot( aes(x=price)) + geom_histogram(color="black", fill="yellow", lwd=0.3)+
   ggtitle("Histograma del precio") +ylab("Frecuencia")

Observamos que la distribución de la variable precio presenta un alto sesgo positivo, por lo que le vamos a aplicar una transformación logarítmica para ver si así consegimos una distribución más simétrica.

p1=d_training %>% ggplot( aes(x= price)) + geom_histogram(aes(y=..density..), bins=30, colour="black", fill="yellow") +
  geom_density(alpha=.5, fill="blue") + ggtitle("Precio")
p2=d_training %>% mutate(log10_p=log(price)) %>%  ggplot( aes(x=log10_p)) + geom_histogram(aes(y=..density..), bins=30, colour="black", fill="yellow") +
  geom_density(alpha=.5, fill="blue") +  ggtitle("Logaritmo del precio")
grid.arrange(p1,p2,nrow=2)

Ahora obtenemos una distribución que a primera vista es como una campana de gauss, de modo que podríamos afirmar que dicha variable sigue una distribución normal. Por ello, a partir de ahora trabajaremos con la variable log(price) en vez de con price.

d_training$price=log(d_training$price)

Bedrooms

d1=d_training %>%  ggplot( aes(x=bedrooms)) + geom_histogram(colour="black", bins =30,fill="tomato")+
  ylab("Frecuencias")

d2=d_training %>%  ggplot( aes(x=as.factor(bedrooms), y=price, fill=as.factor(bedrooms))) + geom_boxplot()+
  labs(x="bedrooms")+theme(legend.position="none")

d3=ggplot(d_training, aes(x = "", y = bedrooms)) + stat_boxplot(geom = "errorbar",  width = 0.3) +
  geom_boxplot(fill = "orange", outlier.colour = "red", alpha = 0.9) +
  xlab("")
grid.arrange(d1,d2,d3,nrow=2) 

En los gráficos representados podemos observar que según aumenta el número de habitaciones aumenta el precio de las casas. También podemos ver posibles outliers, como las casas con 0 y 6 o más habitaciones. Estudiemoslas.

sum(d_training$bedrooms == 0)
## [1] 13

Vemos que hay 13 casas con cero habitaciones. Aunque no es frecuente hay casas, llamadas estudios, en las cuales todo está en una misma estancia. Por ello no consideramos estas observaciones como datos erróneos y las seguimos manteniendo en nuestro dataset.

sum(d_training$bedrooms >=6) 
## [1] 230

Hay 228 casas con 6 y más habitaciones. Al ser un número elevado podemos intuir que tener alto número de habitaciones es algo normal por la zona.

sum(d_training$bedrooms > 8)
## [1] 8

Sin embargo, solo hay 7 casas con 9 o más habitaciones. Esto ya es algo más peculiar, de modo que necesitarán un estudio más profundo.

En principio comprobamos en el mapa la hubicación de estas casas, por si se diese el caso de que están cerca. Si esto fuese así podríamos pensar que es propio de la zona.

Observamos que las casas no están en un misma zona. Lo siguiente que haremos será ver el número de habitaciones que tienen las casas de alrededor.

Comprobando el número de habitaciones de las casas cercanas a nuestros posibles outliers observamos que estas tienen entre 2 y 3 habitaciones. Por ello, el que haya casas por esas zonas con altos números de dormitorios es un poco extraño. Ya que el número de observaciones que tienen más de 8 habitaciones respecto al total de los datos de entrenamiento es muy poco representativo las eliminaremos del estudio.

d_training = d_training[d_training$bedrooms < 8,]

Bathrooms

Observemos los posibles valores que toma esta variable:

table(d_training$bathrooms)
## 
##    0  0.5 0.75    1 1.25  1.5 1.75    2 2.25  2.5 2.75    3 3.25  3.5 3.75    4 
##   10    4   47 2667    5 1007 2135 1341 1417 3779  807  519  410  514  111   90 
## 4.25  4.5 4.75    5 5.25  5.5 5.75    6 6.25  6.5 6.75 7.75 
##   56   69   19   14    9    5    3    3    1    2    1    1

La variable Bathrooms puede tomar valores decimales de 0.25 en 0.25. Y va desde 0 a 8. Esto es debido a que el número de baños se contabiliza por las piezas y cada baño completo tiene 4 piezas.

  • Baño (4 piezas) -> Inodoro, lavabo, bañera y ducha.(1 unidad) - Baño completo
  • Baño (3 piezas) -> Inodoro, lavabo y ducha. (0.75 unidad) - Baño con ducha
  • Baño (2 piezas) -> Inodoro y lavabo. (0.5 unidad) - Aseo
  • Baño (1 pieza) -> Inodoro (0.25 unidad) - Aseo

Vamos a agrupar los baños del siguiente modo:

  • Valores: 0.75, 1 == 1 baño
  • Valores: 1.25, 1.5, 1.75, 2 == 2 baños
  • Valores: 2.25, 2.5, 2.75, 3 == 3 baños

Y así hasta tener 8 baños.

Sin embargo, las casas cuya variable baño tome un valor menor de 0.75 no las consideraremos, ya que una casa que no tenga como mínimo un ducha o un lavabo no es algo lógico. Dependiendo del número de casos en los que ocurra esto, los tomaremos como datos faltantes o simplemente errores de mediciones y si no son muchas las elimiaremos del estudio.

sum(d_training$bathrooms<0.75)
## [1] 14

Obtenemos 14 casas. Que es una número irrelevante frente al total de casos que tenemos. Por lo que eliminaremos estos casos del estudio.

d_training=d_training[-which(d_training$bathrooms<0.75),]

A continuación lo que haremos será cambiar la variable para que considere número de baños en vez de estancias.

d_training$bathrooms=ceiling(d_training$bathrooms)
ggplot(d_training, aes(x = "", y = bathrooms)) +
  stat_boxplot(geom = "errorbar", width = 0.3) +
  geom_boxplot(fill = "orange", outlier.colour = "red", alpha = 0.9) + 
  ggtitle("Boxplot bathrooms") + 
  xlab("") +   
  coord_flip() 

table(d_training$bathrooms) 
## 
##    1    2    3    4    5    6    7    8 
## 2714 4488 6522 1125  158   20    4    1

Podemos observar que hay dos casas con 8 baños y 4 casas con 7 baños. Debido a que esto no es algo común las analizaremos.´

Estudiamos las casas con 7 baños.

d_training[d_training$bathrooms==7,] 
##               id       date    price bedrooms bathrooms sqft_living sqft_lot
## 5667  1924059029 2014-06-17 15.35624        5         7        9640    13068
## 10182 2303900035 2014-06-11 14.87607        5         7        8670    64033
## 14398  424069279 2015-03-28 13.98102        6         7        6260    10955
## 15049 2524069097 2014-05-09 14.62149        5         7        7270   130017
##       floors waterfront view condition grade sqft_above sqft_basement yr_built
## 5667       1          1    4         3    12       4820          4820     1983
## 10182      2          0    4         3    13       6120          2550     1965
## 14398      2          0    0         3    11       4840          1420     2007
## 15049      2          0    0         3    12       6420           850     2010
##       yr_renovated zipcode     lat     long sqft_living15 sqft_lot15
## 5667          2009   98040 47.5570 -122.210          3270      10454
## 10182         2003   98177 47.7295 -122.372          4140      81021
## 14398            0   98075 47.5947 -122.039          2710      12550
## 15049            0   98027 47.5371 -121.982          1800      44890

Nos concuerda el número de baños con los pies cuadrados de la casa, por lo que las conservaremos en nuestro conjunto de datos.
Ahora estudiaremos la casa con 8 baños:

d_training[d_training$bathrooms==8,]
##              id       date    price bedrooms bathrooms sqft_living sqft_lot
## 6467 9208900037 2014-09-19 15.74486        6         8        9890    31374
##      floors waterfront view condition grade sqft_above sqft_basement yr_built
## 6467      2          0    4         3    13       8860          1030     2001
##      yr_renovated zipcode     lat    long sqft_living15 sqft_lot15
## 6467            0   98039 47.6305 -122.24          4540      42730

Observamos un caso que tiene más baños que habitaciones y otro caso que tiene demasiados baños y habitaciones respecto a los pies cuadrados de la casa. Lo cual no tiene mucho sentido. Por lo tanto, eliminamos estas dos observaciones del estudio.

d_training=d_training[-which(d_training$bathrooms==8),]

Agrupamos las casas según el número de baños y realizamos un boxplot para ver si hay diferencias o no en el precio.

d_training %>%  ggplot(aes(x=as.factor(d_training$bathrooms), y=price, fill=as.factor(d_training$bathrooms))) + geom_boxplot()+
  labs(x="bathrooms_group")+theme(legend.position="none")

Se puede ver claramente como según van aumentando el número de baños, aumenta el precio de la casa.

Sqft_living

Veamos su comportamiento frente al precio:

l1<-d_training %>% ggplot(aes(x=sqft_living)) + 
  geom_histogram(aes(y=..density..), bins=30, colour="black", fill="yellow") + 
  geom_density(alpha=.3, fill="blue")

l2<-d_training %>% ggplot(aes(sqft_living, price)) +
  geom_point(alpha = 0.5,col="blue") +
  geom_smooth(se = F, method = "lm", color = "red") +
  scale_y_continuous(breaks = seq(0,8000000, by = 1000000)) 

grid.arrange(l1,l2, nrow=1)

Viendo el scaterplot de los pies cuadrados de la casa frente al precio (gráfico derecho) podemos ver que hay una relación creciente, a medida que aumentan los pies cuadrados de las casas aumenta el precio. Además, en el gráfico izquierdo podemos observar como la distribución de la variable presenta un alto sesgo positivo. Debido a que sería conveniente transformarla para que la distribución de valores fuese más homogénea, le aplicamos una transformación logarítmica.

d_training$sqft_living= log(d_training$sqft_living)

l1log <- ggplot(d_training, aes(x=sqft_living)) + geom_histogram(aes(y=..density..), bins=30, colour="black", fill="white") + geom_density(alpha=.3, fill="#E1AF00")

l2log <- ggplot(d_training, aes(sqft_living, price)) +
  geom_point(alpha = 0.5,col="blue") +
  geom_smooth(se = F, method = "lm", color = "red") +
  scale_y_continuous(breaks = seq(0,8000000, by = 1000000))

grid.arrange(l1log,l2log, nrow=1)

Se puede observar cómo ahora los datos se han normalizado (la distribución presenta la forma de una campana de Gauss).

Basement y above

A continuación vamos a comprobar que la variable above más la variable basement coincide con la variable living. A partir de aquí, vamos a deducir que casas tienen sótano y si el tenerlo o no afecta en el precio.

Total=log(d_training$sqft_above+d_training$sqft_basement)
all(Total==d_training$sqft_living) 
## [1] TRUE
#comprobamos que efectivamente sqft_above+sqft_basement=sqft_living
a=cbind(Total,living=d_training$sqft_living,liv15=log(d_training$sqft_living15),lot=log(d_training$sqft_lot),lot15=(d_training$sqft_lot15))
head(a)
##         Total   living    liv15       lot  lot15
## [1,] 7.073270 7.073270 7.200425  8.639411   5650
## [2,] 7.851661 7.851661 7.432484  8.887653   7639
## [3,] 7.580700 7.580700 7.215240  8.517193   5000
## [4,] 7.426549 7.426549 7.495542  8.997147   7503
## [5,] 8.597851 8.597851 8.468003 11.532042 101930
## [6,] 7.544332 7.544332 7.779049  8.788746   7570
d_training$basement=ifelse(d_training$sqft_living-log(d_training$sqft_above)==0,"0","1") #0 es que no tiene sotano
d_training$basement=as.factor(d_training$basement)

ggplot(d_training, aes(x=basement, y=price, fill=basement)) + geom_boxplot()+
  labs(x="Tener sótano (1) o no tenerlo (0)")+
  scale_y_continuous(labels = scales::dollar,n.breaks = 15)

Se ve un pequeño incremento en el precio en las casas que tienen sótano. Pero no es perceptible, por lo que no consideraremos el tener o no sótano relevante en el modelo.

Living15

Para comprobar si podemos explicar la variable sqft_living15 a partir de sqft_living realizaremos una recta de regresión:

head(a)
##         Total   living    liv15       lot  lot15
## [1,] 7.073270 7.073270 7.200425  8.639411   5650
## [2,] 7.851661 7.851661 7.432484  8.887653   7639
## [3,] 7.580700 7.580700 7.215240  8.517193   5000
## [4,] 7.426549 7.426549 7.495542  8.997147   7503
## [5,] 8.597851 8.597851 8.468003 11.532042 101930
## [6,] 7.544332 7.544332 7.779049  8.788746   7570
cbind(mean(a[,3]-a[,2]),mean(a[,5]-a[,4]))
##              [,1]     [,2]
## [1,] -0.009071659 12714.12
summary(lm(log(d_training$sqft_living15)~d_training$sqft_living, data=d_training))
## 
## Call:
## lm(formula = log(d_training$sqft_living15) ~ d_training$sqft_living, 
##     data = d_training)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.65892 -0.13088  0.00365  0.14332  1.13938 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            3.211374   0.031654   101.5   <2e-16 ***
## d_training$sqft_living 0.573466   0.004186   137.0   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2183 on 15029 degrees of freedom
## Multiple R-squared:  0.5553, Adjusted R-squared:  0.5553 
## F-statistic: 1.877e+04 on 1 and 15029 DF,  p-value: < 2.2e-16

Realizamos un contraste de hipótesis en el cual la hipótesis nula es que la variable no es relevante para el modelo. Esto es, no podemos explicar la variable respuesta (sqft_living15) en función de la explicativa (sqft_living). Al obtener un p.valor menor que 0.05 rechazamos H0 y por consiguiente podemos explicar living15 a partir de living. (Aún así no obtenemos un R^2 ajustado bueno que digamos).

Veamos que variable de entre estas dos explica mejor el precio:

summary(lm(d_training$price~d_training$sqft_living15, data=d_training)) 
## 
## Call:
## lm(formula = d_training$price ~ d_training$sqft_living15, data = d_training)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.88918 -0.30016 -0.00795  0.26023  1.81985 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              1.211e+01  1.037e-02 1167.04   <2e-16 ***
## d_training$sqft_living15 4.736e-04  4.928e-06   96.11   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4145 on 15029 degrees of freedom
## Multiple R-squared:  0.3806, Adjusted R-squared:  0.3806 
## F-statistic:  9237 on 1 and 15029 DF,  p-value: < 2.2e-16
summary(lm(d_training$price~d_training$sqft_living, data=d_training)) 
## 
## Call:
## lm(formula = d_training$price ~ d_training$sqft_living, data = d_training)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.10747 -0.29387  0.01419  0.25964  1.31752 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            6.749564   0.056446   119.6   <2e-16 ***
## d_training$sqft_living 0.834292   0.007464   111.8   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3892 on 15029 degrees of freedom
## Multiple R-squared:  0.4539, Adjusted R-squared:  0.4539 
## F-statistic: 1.249e+04 on 1 and 15029 DF,  p-value: < 2.2e-16

Prediciendo el precio según living obtenemos un R² un poco mayor. Luego entre living y living15 nos quedaremos con living.

Sqft_lot

sl0=ggplot(d_training, aes(x = "", y = sqft_lot)) +
  stat_boxplot(geom = "errorbar",  width = 0.3) +
  geom_boxplot(fill = "orange", outlier.colour = "red",alpha = 0.9) + 
   xlab("") +   coord_flip() 

sl1=ggplot(d_training, aes(x=sqft_lot)) + geom_histogram(aes(y=..density..), bins=30, colour="black", fill="tomato") 

sl2=ggplot(d_training, aes(sqft_lot, price)) +
  geom_point(alpha = 0.5) +
  geom_smooth(se = F, method = "lm", color = "blue") +
  scale_y_continuous(breaks = seq(0,8000000, by = 1000000)) 

grid.arrange(sl0,sl1,sl2, nrow=2)

Esta variable no parece que tenga una relación lineal muy clara con price, por lo que no vemos necesario hacer ninguna transformación ya que no la usaremos para la implementación del modelo. Haciendo el estudio de manera análoga para lot15 llegamos a la misma conclusión.

Floors

table(d_training$floors)
## 
##    1  1.5    2  2.5    3  3.5 
## 7397 1328 5755  115  434    2

No se puede tener un número de pisos que no sea entero. Debido a que no tenemos muy claro el significado de estos decimales, supondremos que indica que la planta está inacabada o que hay espacio en la misma que no se está usando. Atendiendo a esta consideración sobre el número de plantas truncaremos la variable.

d_training$floors=trunc(d_training$floors)

De esta forma nos quedamos con 3 posibles pisos. Veamos si esto influye o no en el precio de las casas.

f1=ggplot(d_training, aes(x=as.factor(floors), y=price, fill=as.factor(floors))) + geom_boxplot()+
  labs(x="Waterfront")+
  scale_y_continuous(labels = scales::dollar,n.breaks = 15)
f2=ggplot(d_training, aes(x=as.factor(floors))) + geom_histogram(colour="black", stat ="count",fill="yellow")+
  labs(x="Condición",y="Frecuencia")
grid.arrange(f2,f1, nrow=1)

Vemos que el numero de plantas no influye en el precio de las casas.

Waterfront

La codificación respecto si tiene o no vistas al mar es 0 y 1, respectivamente.

ggplot(d_training, aes(x=waterfront, y=price, fill=waterfront)) + geom_boxplot()+
  labs(x="Waterfront")+
  scale_y_continuous(labels = scales::dollar,n.breaks = 15)

Observamos un incremento notorio en el precio de las casas que tienen vistas al mar.

Yr_renovated

Vamos a hacer un cambio a esta variable. Y la cambiaremos a sí o no en función de si ha sido o no renovada. Codificamos 0 como no renovada y 1 como renovada, y comprobamos si el que la casa haya sido o no renovada influye en el precio.

d_training$yr_renovated=ifelse(d_training$yr_renovated==0,"0", "1") 
d_training$yr_renovated=as.factor(d_training$yr_renovated)

ggplot(d_training, aes(x=yr_renovated, y=price, fill=yr_renovated)) + geom_boxplot()+
  labs(x="No renovada VS Renovada")+theme(legend.position="none")+
  scale_y_continuous(labels = scales::dollar,n.breaks = 15)

Veamos que sí influye en el precio el que la casa haya sido o no renovada.

Yr_built

Las casas se construyeron entre 1900 y 2015. Veamos si influye esto en el precio

summarise(d_training, añoMin=min(d_training$yr_built),añoMax= max(d_training$yr_built))
##   añoMin añoMax
## 1   1900   2015
ggplot(d_training, aes(x=yr_built)) + geom_histogram(colour="black",bins=30, fill="#E1AF00")+
  labs(title="Histograma del precio y año de construcción", x="Año construida",y="Precio")

Vemos un claro aumento en el precio de las casas construidas en los últimos años.

Condition

table(d_training$condition)
## 
##    1    2    3    4    5 
##   23  119 9767 3951 1171
c1=ggplot(d_training, aes(x=condition, y=price, fill=condition)) + geom_boxplot()+
  ylab("Precio")+theme(legend.position="none")+
  scale_y_continuous(labels = scales::dollar,n.breaks = 15)

c2=ggplot(d_training, aes(x=condition)) + geom_histogram(colour="black", stat ="count",fill="yellow")+
  ylab("Frecuencia") 

grid.arrange(c1,c2, nrow=1)

Observamos que sí es una variable significativa.

View

v1=ggplot(d_training, aes(x=view, y=price, fill=view)) + geom_boxplot()+
  ylab("Precio")+
  scale_y_continuous(labels = scales::dollar,n.breaks = 15) 

v2=ggplot(d_training, aes(x=view)) + geom_histogram(colour="black", stat ="count",fill="yellow")+
 ylab("Frecuencia")
grid.arrange(v1,v2, nrow=1)

Observamos que si la casa tiene mejores vistas su precio es mayor. Sin embargo, no hay muchas casas con buenas vistas.

Grade

La variable grade mide el nivel de construcción y diseño de la casa en un rango del 1 al 13. Entre los valores 1 y 3 indica que no llega a la construcción y el diseño de edificios, 7 que tiene un nivel promedio, y entre 11 y 13 indica que tiene un alto nivel de calidad de construcción y diseño.

En función de la descripción de la variable tomada, vamos a agruparla en 4 para reducir el número de sus posibles niveles.

  • Inacabada = 0 a 3
  • Aceptable = 4 a 7
  • Buena = 8 a 10
  • Excelente = 11 a 13
d_training$grade <- cut(d_training$grade, breaks = c(0,3,7,10,13), labels = c("Inacabada","Aceptable","Buena","Excelente"))

ggplot(d_training, aes(x=grade, y=price, fill=grade)) + geom_boxplot()+
  labs(x="Grado de construcción",y="Precio")+
  scale_y_continuous(labels = scales::dollar,n.breaks = 15)

El grado de construcción sí que afecta al precio, y bastante como podemos observar.

Zipcode

Realizamos un boxplot diferenciando con colores los distintos valores de zipcode:

ggplot(d_training, aes(x=zipcode, y=price, fill=zipcode)) + geom_boxplot()+
  ylab("Precio")+ guides(fill=FALSE) +
  scale_y_continuous(labels = scales::dollar,n.breaks = 15) +
  theme(
        axis.text.x = element_text(angle=90,size=7)
     )

A simple vista no se observa ningún patrón de precios dependiendo del área donde esté ubicada la casa.

6. Ajuste, interpretación y diagnosis del modelo de regresión lineal

Nuestro objetivo es predecir el precio de las casas en el condado de King y comprender qué factores son los responsables de un valor de propiedad más alto.

Antes de nada eliminaremos del conjunto de dato la variable “identificador” ya que no es relevante para predecir el precio de la vivienda. No es más que un identificador de la casa, como sus nombre por así decirlo .

d_training=select(d_training,-id)

A continuación ajustaremos una recta de regresión para predecir la variable respuesta precio en función de las variables explicativas del modelo.

Veamos en principio la colinealidad entre las varibles cuantitativas

d_train_cuant=select(d_training,-date,-waterfront,-view,-condition,-basement,-grade,-yr_built,-yr_renovated,-zipcode) #seleccionamos las variables cuantitativas

correlaciones= cor(d_train_cuant) #correlación entre las variables cuantitativas
correlaciones
##                    price     bedrooms  bathrooms sqft_living     sqft_lot
## price         1.00000000  0.355195293 0.52864672  0.67374379  0.090451314
## bedrooms      0.35519529  1.000000000 0.50620174  0.64462664  0.027010099
## bathrooms     0.52864672  0.506201745 1.00000000  0.74934905  0.067675948
## sqft_living   0.67374379  0.644626641 0.74934905  1.00000000  0.134933646
## sqft_lot      0.09045131  0.027010099 0.06767595  0.13493365  1.000000000
## floors        0.28801889  0.163688858 0.55381091  0.36347719 -0.009682656
## sqft_above    0.60316447  0.499271684 0.66635231  0.83440895  0.171245781
## sqft_basement 0.31084025  0.301591608 0.25086854  0.41639899 -0.001037189
## lat           0.45036426 -0.008473969 0.01652289  0.03893912 -0.094040091
## long          0.04919025  0.144571652 0.24061825  0.25933751  0.227252958
## sqft_living15 0.61696436  0.410253085 0.56037916  0.73464487  0.135475873
## sqft_lot15    0.08100238  0.026356736 0.06888082  0.15188713  0.696687244
##                     floors   sqft_above sqft_basement          lat        long
## price          0.288018888  0.603164473   0.310840247  0.450364264  0.04919025
## bedrooms       0.163688858  0.499271684   0.301591608 -0.008473969  0.14457165
## bathrooms      0.553810915  0.666352306   0.250868537  0.016522892  0.24061825
## sqft_living    0.363477193  0.834408954   0.416398988  0.038939116  0.25933751
## sqft_lot      -0.009682656  0.171245781  -0.001037189 -0.094040091  0.22725296
## floors         1.000000000  0.518092745  -0.242020536  0.032012208  0.15984579
## sqft_above     0.518092745  1.000000000  -0.062427518 -0.002549824  0.34877063
## sqft_basement -0.242020536 -0.062427518   1.000000000  0.112435799 -0.15086633
## lat            0.032012208 -0.002549824   0.112435799  1.000000000 -0.14519166
## long           0.159845794  0.348770625  -0.150866335 -0.145191661  1.00000000
## sqft_living15  0.291342424  0.728462262   0.197486776  0.045002865  0.33933155
## sqft_lot15    -0.020755199  0.179429883   0.012308039 -0.101471044  0.25137265
##               sqft_living15  sqft_lot15
## price            0.61696436  0.08100238
## bedrooms         0.41025309  0.02635674
## bathrooms        0.56037916  0.06888082
## sqft_living      0.73464487  0.15188713
## sqft_lot         0.13547587  0.69668724
## floors           0.29134242 -0.02075520
## sqft_above       0.72846226  0.17942988
## sqft_basement    0.19748678  0.01230804
## lat              0.04500286 -0.10147104
## long             0.33933155  0.25137265
## sqft_living15    1.00000000  0.18380430
## sqft_lot15       0.18380430  1.00000000
corrplot(correlaciones,type="upper")

Correlaciones con la variable respuesta:

correlaciones[,1] 
##         price      bedrooms     bathrooms   sqft_living      sqft_lot 
##    1.00000000    0.35519529    0.52864672    0.67374379    0.09045131 
##        floors    sqft_above sqft_basement           lat          long 
##    0.28801889    0.60316447    0.31084025    0.45036426    0.04919025 
## sqft_living15    sqft_lot15 
##    0.61696436    0.08100238

Con respecto al precio vemos alta correlación con el sqft_living y sqft_above. De modo que en el modelo futuro se espera que estas variables tengan mayor peso.

det(correlaciones)
## [1] 0.0004694057

El determinante de la matriz de correlaciones entre las variables explicativas es cero. Esto implica que si consideramos todas estas variables tendremos un problema de multicolinealidad. Por ello, para construir la recta de regresión realizaremos una selección de variables a mano.

Según el estudio visual de las variables más la matriz de correlaciones tomamos en principio las siguiente selección de variables:

model1 <- lm(price~bedrooms+bathrooms+sqft_living+waterfront+view+condition+
               grade+yr_renovated+yr_built+lat, data=d_training)
summary(model1)
## 
## Call:
## lm(formula = price ~ bedrooms + bathrooms + sqft_living + waterfront + 
##     view + condition + grade + yr_renovated + yr_built + lat, 
##     data = d_training)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.260 -0.175  0.000  0.166  1.228 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -5.525e+01  9.080e-01 -60.854  < 2e-16 ***
## bedrooms       -3.877e-02  3.387e-03 -11.448  < 2e-16 ***
## bathrooms       6.125e-02  4.250e-03  14.412  < 2e-16 ***
## sqft_living     5.826e-01  9.878e-03  58.978  < 2e-16 ***
## waterfront1     3.549e-01  3.237e-02  10.963  < 2e-16 ***
## view1           1.976e-01  1.783e-02  11.084  < 2e-16 ***
## view2           1.460e-01  1.093e-02  13.362  < 2e-16 ***
## view3           2.308e-01  1.501e-02  15.373  < 2e-16 ***
## view4           3.040e-01  2.278e-02  13.346  < 2e-16 ***
## condition2      8.360e-02  6.220e-02   1.344 0.178935    
## condition3      1.915e-01  5.716e-02   3.350 0.000811 ***
## condition4      2.388e-01  5.715e-02   4.179 2.94e-05 ***
## condition5      2.914e-01  5.752e-02   5.066 4.10e-07 ***
## gradeAceptable -6.251e-01  2.732e-01  -2.288 0.022156 *  
## gradeBuena     -3.898e-01  2.733e-01  -1.426 0.153851    
## gradeExcelente  2.333e-02  2.739e-01   0.085 0.932135    
## yr_renovated1   6.725e-02  1.187e-02   5.667 1.48e-08 ***
## yr_built       -2.630e-03  1.094e-04 -24.041  < 2e-16 ***
## lat             1.458e+00  1.668e-02  87.397  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2729 on 15012 degrees of freedom
## Multiple R-squared:  0.7319, Adjusted R-squared:  0.7316 
## F-statistic:  2277 on 18 and 15012 DF,  p-value: < 2.2e-16

Obtenemos unos coeficientes lógicos y un R² ajustado de 0.7316. En cambio indica que a mayor número de habitaciones menor precio, lo cual es un poco raro.

En general es buen modelo, aunque la variable grade no parece influir mucho. Vemos que pasa si no la introducimos:

model1b=lm(price~bedrooms+bathrooms+sqft_living+waterfront+view+condition+
    yr_renovated+yr_built+lat, data=d_training) #model1 sin el grade
summary(model1b) 
## 
## Call:
## lm(formula = price ~ bedrooms + bathrooms + sqft_living + waterfront + 
##     view + condition + yr_renovated + yr_built + lat, data = d_training)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.23325 -0.18387 -0.00279  0.17633  1.25404 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6.263e+01  9.153e-01 -68.426  < 2e-16 ***
## bedrooms      -6.310e-02  3.571e-03 -17.672  < 2e-16 ***
## bathrooms      9.114e-02  4.483e-03  20.329  < 2e-16 ***
## sqft_living    7.487e-01  9.794e-03  76.440  < 2e-16 ***
## waterfront1    3.585e-01  3.454e-02  10.379  < 2e-16 ***
## view1          2.222e-01  1.902e-02  11.679  < 2e-16 ***
## view2          1.874e-01  1.163e-02  16.117  < 2e-16 ***
## view3          2.874e-01  1.597e-02  17.997  < 2e-16 ***
## view4          3.727e-01  2.425e-02  15.366  < 2e-16 ***
## condition2     6.170e-02  6.640e-02   0.929 0.352731    
## condition3     1.766e-01  6.102e-02   2.894 0.003809 ** 
## condition4     2.147e-01  6.101e-02   3.520 0.000433 ***
## condition5     2.601e-01  6.140e-02   4.235 2.30e-05 ***
## yr_renovated1  6.729e-02  1.265e-02   5.318 1.06e-07 ***
## yr_built      -1.907e-03  1.145e-04 -16.658  < 2e-16 ***
## lat            1.547e+00  1.766e-02  87.572  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2913 on 15015 degrees of freedom
## Multiple R-squared:  0.6944, Adjusted R-squared:  0.6941 
## F-statistic:  2275 on 15 and 15015 DF,  p-value: < 2.2e-16
model1c <- lm(price~bathrooms+sqft_living+waterfront+view+condition+
               grade+yr_renovated+yr_built+lat, data=d_training) #modelo1 sin bedroom
summary(model1c)
## 
## Call:
## lm(formula = price ~ bathrooms + sqft_living + waterfront + view + 
##     condition + grade + yr_renovated + yr_built + lat, data = d_training)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.26112 -0.17537  0.00039  0.16626  1.23274 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -5.563e+01  9.113e-01 -61.042  < 2e-16 ***
## bathrooms       5.489e-02  4.232e-03  12.971  < 2e-16 ***
## sqft_living     5.272e-01  8.648e-03  60.963  < 2e-16 ***
## waterfront1     3.681e-01  3.249e-02  11.330  < 2e-16 ***
## view1           2.051e-01  1.789e-02  11.463  < 2e-16 ***
## view2           1.529e-01  1.096e-02  13.949  < 2e-16 ***
## view3           2.418e-01  1.505e-02  16.071  < 2e-16 ***
## view4           3.141e-01  2.286e-02  13.742  < 2e-16 ***
## condition2      8.534e-02  6.246e-02   1.366 0.171910    
## condition3      1.903e-01  5.741e-02   3.315 0.000919 ***
## condition4      2.364e-01  5.739e-02   4.118 3.84e-05 ***
## condition5      2.887e-01  5.777e-02   4.998 5.86e-07 ***
## gradeAceptable -6.527e-01  2.744e-01  -2.379 0.017370 *  
## gradeBuena     -4.079e-01  2.745e-01  -1.486 0.137278    
## gradeExcelente  1.995e-02  2.751e-01   0.073 0.942200    
## yr_renovated1   7.311e-02  1.191e-02   6.140 8.43e-10 ***
## yr_built       -2.503e-03  1.093e-04 -22.899  < 2e-16 ***
## lat             1.467e+00  1.673e-02  87.695  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2741 on 15013 degrees of freedom
## Multiple R-squared:  0.7296, Adjusted R-squared:  0.7293 
## F-statistic:  2383 on 17 and 15013 DF,  p-value: < 2.2e-16

El modelo llamado “model1b” indica que mientras más habitaciones más baratas las casas… sospechoso y un R² = 0.6941. EL model1c presenta en general unos coeficientes lógicos, aunque los de condi2 y vieew1 nos chocan un poco y un R²=0.7293.

model1cc <- lm(price~bathrooms+sqft_living+waterfront+view+condition+
                yr_renovated+yr_built+lat, data=d_training) # modelo1c sin grade
summary(model1cc)
## 
## Call:
## lm(formula = price ~ bathrooms + sqft_living + waterfront + view + 
##     condition + yr_renovated + yr_built + lat, data = d_training)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.23353 -0.18880 -0.00439  0.18000  1.26612 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6.376e+01  9.225e-01 -69.111  < 2e-16 ***
## bathrooms      8.249e-02  4.503e-03  18.320  < 2e-16 ***
## sqft_living    6.671e-01  8.727e-03  76.438  < 2e-16 ***
## waterfront1    3.806e-01  3.487e-02  10.915  < 2e-16 ***
## view1          2.363e-01  1.920e-02  12.306  < 2e-16 ***
## view2          2.015e-01  1.172e-02  17.195  < 2e-16 ***
## view3          3.094e-01  1.608e-02  19.236  < 2e-16 ***
## view4          3.940e-01  2.447e-02  16.098  < 2e-16 ***
## condition2     6.313e-02  6.708e-02   0.941 0.346670    
## condition3     1.737e-01  6.165e-02   2.817 0.004847 ** 
## condition4     2.091e-01  6.164e-02   3.392 0.000695 ***
## condition5     2.536e-01  6.204e-02   4.087 4.39e-05 ***
## yr_renovated1  7.719e-02  1.277e-02   6.045 1.53e-09 ***
## yr_built      -1.645e-03  1.147e-04 -14.344  < 2e-16 ***
## lat            1.568e+00  1.780e-02  88.111  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2943 on 15016 degrees of freedom
## Multiple R-squared:  0.6881, Adjusted R-squared:  0.6878 
## F-statistic:  2366 on 14 and 15016 DF,  p-value: < 2.2e-16

Obtenemos un R² de 0.6878, lo cual empeoramos. Rechazamos esta selección de variables.

model1a=(lm(price~sqft_living+waterfront+view+condition+
             +grade+yr_renovated+yr_built+lat, data=d_training)) #model1 quitando bathrooms y bedrooms
summary(model1a)#coeficientes lógicos, aunque cond2 y view1 no del todo y R²= 0.7263
## 
## Call:
## lm(formula = price ~ sqft_living + waterfront + view + condition + 
##     +grade + yr_renovated + yr_built + lat, data = d_training)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.20412 -0.17392  0.00029  0.16770  1.22714 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -5.767e+01  9.026e-01 -63.896  < 2e-16 ***
## sqft_living     5.926e-01  7.061e-03  83.926  < 2e-16 ***
## waterfront1     3.720e-01  3.267e-02  11.386  < 2e-16 ***
## view1           2.065e-01  1.799e-02  11.478  < 2e-16 ***
## view2           1.555e-01  1.102e-02  14.111  < 2e-16 ***
## view3           2.467e-01  1.513e-02  16.312  < 2e-16 ***
## view4           3.137e-01  2.299e-02  13.646  < 2e-16 ***
## condition2      8.442e-02  6.281e-02   1.344  0.17898    
## condition3      1.890e-01  5.772e-02   3.275  0.00106 ** 
## condition4      2.348e-01  5.771e-02   4.068 4.77e-05 ***
## condition5      2.961e-01  5.809e-02   5.098 3.48e-07 ***
## gradeAceptable -6.899e-01  2.759e-01  -2.501  0.01240 *  
## gradeBuena     -4.360e-01  2.760e-01  -1.580  0.11419    
## gradeExcelente  6.933e-03  2.766e-01   0.025  0.98001    
## yr_renovated1   9.004e-02  1.190e-02   7.566 4.07e-14 ***
## yr_built       -1.922e-03  1.002e-04 -19.170  < 2e-16 ***
## lat             1.479e+00  1.680e-02  88.057  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2756 on 15014 degrees of freedom
## Multiple R-squared:  0.7266, Adjusted R-squared:  0.7263 
## F-statistic:  2493 on 16 and 15014 DF,  p-value: < 2.2e-16
model2 =lm(price~bathrooms+sqft_living+waterfront+view+condition+
             grade+yr_renovated+lat, data=d_training) #modelo1 sin bedroom ni yr_built
summary(model2)#los coeficientes del grade no son lógicos
## 
## Call:
## lm(formula = price ~ bathrooms + sqft_living + waterfront + view + 
##     condition + grade + yr_renovated + lat, data = d_training)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.22004 -0.17945 -0.00649  0.16520  1.28628 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -64.813514   0.832423 -77.861  < 2e-16 ***
## bathrooms        0.015177   0.003927   3.865 0.000111 ***
## sqft_living      0.546092   0.008757  62.362  < 2e-16 ***
## waterfront1      0.368346   0.033051  11.145  < 2e-16 ***
## view1            0.227822   0.018172  12.537  < 2e-16 ***
## view2            0.181053   0.011078  16.344  < 2e-16 ***
## view3            0.274696   0.015238  18.027  < 2e-16 ***
## view4            0.337649   0.023231  14.534  < 2e-16 ***
## condition2       0.055359   0.063529   0.871 0.383548    
## condition3       0.114851   0.058301   1.970 0.048861 *  
## condition4       0.199371   0.058364   3.416 0.000637 ***
## condition5       0.279342   0.058767   4.753 2.02e-06 ***
## gradeAceptable  -0.722114   0.279092  -2.587 0.009681 ** 
## gradeBuena      -0.506090   0.279220  -1.813 0.069927 .  
## gradeExcelente  -0.071329   0.279833  -0.255 0.798804    
## yr_renovated1    0.158466   0.011504  13.775  < 2e-16 ***
## lat              1.558626   0.016531  94.288  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2788 on 15014 degrees of freedom
## Multiple R-squared:  0.7201, Adjusted R-squared:  0.7198 
## F-statistic:  2415 on 16 and 15014 DF,  p-value: < 2.2e-16
model3= lm(price~bathrooms+sqft_living+waterfront+view+condition+
     grade+yr_built+lat, data=d_training) #modelo1 sin bedroom sin bedroom ni yr_renovated
summary(model3)#coeficientes que no son lógicos y R² = 0.7286 
## 
## Call:
## lm(formula = price ~ bathrooms + sqft_living + waterfront + view + 
##     condition + grade + yr_built + lat, data = d_training)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.20406 -0.17547  0.00073  0.16582  1.22462 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -5.500e+01  9.066e-01 -60.664  < 2e-16 ***
## bathrooms       5.774e-02  4.211e-03  13.710  < 2e-16 ***
## sqft_living     5.281e-01  8.657e-03  61.009  < 2e-16 ***
## waterfront1     3.771e-01  3.250e-02  11.604  < 2e-16 ***
## view1           2.075e-01  1.791e-02  11.585  < 2e-16 ***
## view2           1.528e-01  1.097e-02  13.922  < 2e-16 ***
## view3           2.440e-01  1.506e-02  16.198  < 2e-16 ***
## view4           3.189e-01  2.287e-02  13.943  < 2e-16 ***
## condition2      8.701e-02  6.254e-02   1.391 0.164156    
## condition3      1.973e-01  5.746e-02   3.433 0.000598 ***
## condition4      2.380e-01  5.746e-02   4.141 3.48e-05 ***
## condition5      2.877e-01  5.784e-02   4.974 6.61e-07 ***
## gradeAceptable -6.492e-01  2.747e-01  -2.363 0.018132 *  
## gradeBuena     -4.032e-01  2.748e-01  -1.467 0.142422    
## gradeExcelente  2.096e-02  2.754e-01   0.076 0.939333    
## yr_built       -2.713e-03  1.039e-04 -26.103  < 2e-16 ***
## lat             1.462e+00  1.673e-02  87.393  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2744 on 15014 degrees of freedom
## Multiple R-squared:  0.7289, Adjusted R-squared:  0.7286 
## F-statistic:  2523 on 16 and 15014 DF,  p-value: < 2.2e-16

Teniendo en cuenta el R² y los coeficientes de las variables seleccionadas en cada modelo anterior nos quedaremos con el modelo denominado “model1a” :
model1a=lm(price~sqft_living+waterfront+view+condition+
+grade+yr_renovated+yr_built+lat, data=d_training) y veamos cómo de bien o mal precide el precio de las casas del conjunto train.

Colinealidad en el modelo Para ver si existe colinealidad entre las variables seleccionadas en nuestro modelo seleccionado vamos a usar el Factor de Inflación de la Varianza (VIF),

Los límites de referencia que se suelen emplear son:

  • VIF = 1: Ausencia total de colinealidad.
  • 1 < VIF < 5: La regresión puede verse afectada por cierta colinealidad.
  • 5 < VIF < 10: Causa de preocupación.
vif(model1a)
##                  GVIF Df GVIF^(1/(2*Df))
## sqft_living  1.785524  1        1.336235
## waterfront   1.506865  1        1.227544
## view         1.678938  4        1.066914
## condition    1.282599  4        1.031600
## grade        1.937298  3        1.116519
## yr_renovated 1.123988  1        1.060183
## yr_built     1.734619  1        1.317050
## lat          1.075865  1        1.037239

Observamos que los valores se encuentran entre 1 y 2 y por consiguiente no nos indican la presencia de colinealidad entre variables.

Análisis de residuos

Una vez estimado el modelo de regresión y obtenido los residuos hay que comprobar si las hipótesis que se han utilizado para construirlo se pueden asumir como ciertas o no. Si no lo son, habrá que modificar el modelo para adaptarlo a los datos observados.

Para hacer el modelo se asumen 4 hipótesis: * 1 Linealidad: La relación entre la variable respuesta y las variables explicativas es lineal. * 2 Homocedasticidad: La variabilidad del error es constante. Esto es, el error sigue una distribución normal con varianza desconocida, pero todas iguales y constantes. * 3 Normalidad: Las perturbaciones (el error) siguen una distribución normal. * 4 Independencia: Las perturbaciones son independientes entre sí.

par(mfrow=c(2,2))
plot(model1a)

Linealidad

Se está analizando si existe una relación lineal entre la variable respuesta y las variables explicativas. Esto lo comprobamos en la primera de las 4 gráficas que tenemos.

Si se verifica la hipótesis de linealidad, esta gráfica debería de presentar una simetría respecto al eje horizontal. Como vemos en nuestro gráfico si se verifica la linealidad.

Normalidad

Para comprobar la normalidad de los residuos hay que hacer un gráfico cuantil-cuantil (un Q-Q plot) de los residuos, que se corresponde con la gŕafica de la esquina superior derecha.

A simple vista podŕiamos decir que a pesar de las colas si siguen una distribución normal. Comprobémoslo usando un contraste de normalidad:

residuos=resid(model1a)
ks.test(residuos, 'pnorm')
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  residuos
## D = 0.27828, p-value < 2.2e-16
## alternative hypothesis: two-sided

Se contrasta la hipótesis nula de que los residuos del modelo se distribuyen según una Normal. Al obtener un p.valor menor que 0.05 rechazamos H0 y por consiguiente se tiene que los residuos no siguen una distribución Normal. Sin embargo, debido a que los modelos lineales en general son robustos a la ligera falta de normalidad asumiremos la hipótesis de normalidad como cierta y proseguiremos con el estudio.

Homocedasticidad

Para comprobar la homocedasticidad hay que hacer un gráfico de dispersión cuantil-cuantil de los residuos. Lo ideal es que la recta sea parecida a y=0. Debemos obtener una línea recta horizontal. Si no fuese así lo que nos está mostrando es cierta dependencia entre la magnitud de los errores y la predicción.

El gŕafico inferior derecho es un gráfico de leverage, de apalancamiento, que sirve para comprobar si hay alguna observación que sea demasiado influyente en la construcción de los coeficientes del modelo. Si los puntos están agrupados y sin sobrepasar las curvas de nivel de la distancia de cook no habrá ningún problema. Si hay algún dato muy influyente, el gráfico te da un número indicando cuál es exactamente y tendríamos que estudiarlas para ver si son un error o un dato verdaderamente distinto al resto (un dato que sigue un patrón muy distinto al resto) e investigar su causa. El estudio de estos valores atípicos puede ser interesante porque pueden modificar los coeficientes del modelo

7. Predicción sobre los datos de testing. Evaluación del modelo

Una vez construido el modelo predictor, comprobemos como de bueno o malo o es. Esto lo haremos sobre el conjunto test que tenemos. Vamos a comprobar si con el modelo elegido obtenemos una buena predicción del precio.

Primero, realizaremos las mismas transformaciones, agrupaciones y categorizaciones que realizamos con los datos de train.

d_testing$price <- log(d_testing$price)

d_testing$bathrooms=ceiling(d_testing$bathrooms)

d_testing$sqft_living<- log(d_testing$sqft_living)

d_testing$grade <- cut(d_testing$grade, breaks = c(0,3,7,10,13), labels = c("Inacabada","Aceptable","Buena","Excelente"))


d_testing$basement=ifelse(d_testing$sqft_living-log(d_testing$sqft_above)==0,"0","1") 
d_testing$basement=as.factor(d_testing$basement)

d_testing$yr_renovated=ifelse(d_testing$yr_renovated==0,"0", "1") 
d_testing$yr_renovated=as.factor(d_testing$yr_renovated)

Ahora veamos cómo predice el modelo:

pred=predict(model1a,d_testing[,-3])
New=as.data.frame(cbind(Prediccion=exp(pred),Observ=exp(d_testing[,3]), Dif=exp(pred)-exp(d_testing[,3])))
summary(New)
##    Prediccion          Observ             Dif          
##  Min.   : 147591   Min.   :  78000   Min.   :-3443298  
##  1st Qu.: 337479   1st Qu.: 321625   1st Qu.:  -80093  
##  Median : 450798   Median : 450000   Median :   -1071  
##  Mean   : 515613   Mean   : 538148   Mean   :  -22416  
##  3rd Qu.: 629282   3rd Qu.: 645000   3rd Qu.:   71015  
##  Max.   :4699766   Max.   :7700000   Max.   : 1514141  
##  NA's   :22                          NA's   :22

Obtenemos que el valor absoluto de la media de la diferencia entre el precio predicho y el observado es de 22416. Debido a que tenemos un error pequeño de predicción podemos considerar que tenemos un buen modelo.

plot(exp(d_testing$price),exp(pred),col="blue", xlab = "Precio según el modelo",
     ylab = "Precio observado", main="Representación del precio predicho frente al observado")